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Scientific  Progress 

Fundamental  research  in  the  area  of  solid-state  devices  is  a  corner  stone  and  an  essential  requirement  for  the  development  of 
future  electronic  systems  with  relevance  to  military  defense.  In  order  to  offer  superior  capabilities,  these  systems  will  need  to 
operate  at  much  higher  speeds  and  frequencies,  have  greatly  increased  functionality,  and  have  much  higher  levels  of 
integration 

than  present  day  technology  provides.  Resonant  tunneling  diodes  (RTDs)  have  been  of  interest  to  the  physics  and  engineering 
communities  since  the  1980s  because  they  are  inherently  high-speed  devices.  For  example,  an  RTD-based  fundamental 
oscillator  achieved  the  highest  frequency  generated  by  a  solid-state  source  (>  700  GHz)  in  1991  [1],  The  simulation  a 
resonant-tunneling  diode  (RTD)  can  be  made  in  different  ways.  The  first  successful  approach  was  made  solving  the 
Schrodinger  equation  and  calculating  the  transmission  of  a  diode.  However  there  is  a  problem  how  to  include  the  scattering 
processes  in  quantum  mechanical  simulation  a  resonant-tunneling  diode  .  The  Green  function  method  has  been  one  of 
alternatives.  Another  alternative  was  to  utilize  the  Wigner  function.  The  Wigner  function  approach  to  the  simulation  of  a 
resonant-tunneling  diode  offers  many  advantages.  In  the  limit  of  the  classical  physics  the  Wigner  equation  turns  into  the 
Boltzmann  equation  .  The  scattering  processes  can  be  presented  in  the  similar  way  as  in  the  Boltzmann  equation  .  A  short 
review  of  simulations  made  by  different  authors  using  the  Wigner  formalism  can  be  found  in  [2].  The  Wigner  function  equation 
was  first  employed  in  quantum  device  simulation  by  Frensley  [3].  Later,  Kluksdahl  et  al  [4]  incorporated  Poisson's  equation 
(PE)and  applied  the  model  to  the  study  of  RTD  with  self-consistent  potentials.  They  have  chosen  Lax-Wendroff  explicit  time 
differencing.  Jensen,  and  Buot  [5]  have  used  Cayley  time  evolution  operator. 

Recent  work  has  focused  on  solving  the  time-dependent  Wigner-Poisson  equations  to  look  for  oscillatory  solutions.  Several 
groups  have  studied  dynamical  instabilities  and  undertaken  the  simulation  the  RTD  time-dependent  behavior  using  the  discrete 
Wigner  function  approach  to  look  for  oscillatory  solutions  [5-1 0],  However  there  are  also  difficulties  in  applications  of  the 
Wigner  formalism.  It  is  important  to  recognize  that  the  Wigner  function  is  defined  over  both  real  and  momentum  spaces  (i.e., 
which  differs  from  all  other  quantum  mechanical  formalisms)  and  this  significantly  increases  the  theoretical  difficulty  associated 
with  quantifying  the  exact  requirements  on  the  numerical  discretization  so  as  to  insure  that  all  necessary  physical  laws  are 
satisfied  (e.g.,  charge  continuity,  momentum  balance,  detailed  balance  of  the  equilibrium  state,  and  stability  of  non-equilibrium 
states).  Furthermore,  the  more  refined  numerical  solutions  of  the  Wigner-Poisson  model  yields  results  in  disagreement  with 
prior  results  that  utilized  more  coarse  grained  discretizations. 

This  proposed  research  project  studies  the  discretization  issues  associated  with  the  stability  and  physical  accuracy  of  the 
numerical  solutions  to  the  Wigner-Poisson  modeling  of  resonant  tunneling  diodes.  Here  the  goal  is  to  establish  the  qualitative 
and  quantitative  requirements  for  the  phase  space  discretizations  that  will  be  able  to  generate  physically  accurate  numerical 
solutions  to  the  Wigner-Poisson  model  for  open  systems.  The  goals  of  this  project  is  also  to  study  stability  and  physical 
accuracy  of  numerical  solutions  to  the  Wigner-Poisson  model. 

All  authors  used  the  discrete  Wigner  function  which  is  periodic  in  momentum  space.  The  periodicity  was  proved  from  the 
Fourier  transform  of  the  density  matrix.  The  inverse  Fourier  transform  provides  us  with  the  important  proof  that  the  number  of 
spatial  intervals  is  equal  to  the  number  of  intervals  in  the  momentum  space.  In  addition  we  obtain  that  the  step  size  in  the 
momentum  space  does  not  depend  on  the  number  of  intervals.  As  a  result  the  number  of  relevant  intervals  in  the  momentum 
space  does  not  also  depend  on  the  total  number  of  intervals. 

The  analysis  of  the  stability  shows  that  Greg  Recine  has  used  an  unstable  algorithm. 
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Abstract. 

The  Wigner  formalism  has  previously  been  applied  by  numerous  groups  to  analyze 
the  steady-state  behavior  of  particle  transport  within  open  quantum  systems.  In  particular,  the 
Wigner-Poisson  model  has  been  used  as  a  tool  to  study  electron  transport  through  double -barrier 
resonant  tunneling  diodes  (RTDs)  with  open  boundary  conditions.  The  goals  of  this  project  is  to 
study  stability  and  physical  accuracy  of  numerical  solutions  to  the  Wigner-Poisson  model. 
Numerous  authors  who  simulated  RTD  with  the  Wigner-Poisson  model  assumed  that  the  discrete 
Wigner  function  is  periodic  in  momentum  space.  The  periodicity  follows  from  the  Fourier 
transform  of  the  density  matrix.  The  inverse  Fourier  transform  provides  us  with  the  important 
proof  that  the  number  of  spatial  intervals  is  equal  to  the  number  of  intervals  in  the  momentum 
space.  In  addition  we  obtain  that  the  step  size  in  the  momentum  space  does  not  depend  on  the 
number  of  intervals.  As  a  result  the  number  of  relevant  intervals  in  the  momentum  space  does  not 
also  depend  on  the  total  number  of  intervals.  Stability  of  algorithms  used  by  different  authors  is 
investigated.  The  analysis  of  the  stability  shows  that  Greg  Recine  has  used  unstable  algorithm.  . 

1. Introduction. 


Fundamental  research  in  the  area  of  solid-state  devices  is  a  corner  stone  and  an  essential 
requirement  for  the  development  of  future  electronic  systems  with  relevance  to  military  defense. 
In  order  to  offer  superior  capabilities,  these  systems  will  need  to  operate  at  much  higher  speeds 
and  frequencies,  have  greatly  increased  functionality,  and  have  much  higher  levels  of  integration 
than  present  day  technology  provides.  Resonant  tunneling  diodes  (RTDs)  have  been  of  interest  to 
the  physics  and  engineering  communities  since  the  1980s  because  they  are  inherently  high-speed 
devices.  For  example,  an  RTD-based  fundamental  oscillator  achieved  the  highest  frequency 
generated  by  a  solid-state  source  (>  700  GFlz)  in  1991  [1.1]. 

The  simulation  a  resonant-tunneling  diode  (RTD)  can  be  made  in  different  ways.  The 
first  successful  approach  was  made  solving  the  Schrodinger  equation  and  calculating  the 
transmission  of  a  diode.  However  there  is  a  problem  how  to  include  the  scattering  processes  in 
quantum  mechanical  simulation  a  resonant-tunneling  diode  .  The  Green  function  method  has  been 
one  of  alternatives.  Another  alternative  was  to  utilize  the  Wigner  function  . 

The  Wigner  function  approach  to  the  simulation  of  a  resonant-tunneling  diode  offers 
many  advantages.  In  the  limit  of  the  classical  physics  the  Wigner  equation  turns  into  the 
Boltzmann  equation  .  The  scattering  processes  can  be  presented  in  the  similar  way  as  in  the 
Boltzmann  equation  .  A  short  review  of  simulations  made  by  different  authors  using  the  Wigner 
formalism  can  be  found  in  [2]. 

The  Wigner  function  equation  was  first  employed  in  quantum  device  simulation 


by  Frensley  [3].  Later,  Kluksdahl  et  al  [4]  incorporated  Poisson's  equation  (PE)  and  applied  the 
model  to  the  study  of  RTD  with  self-consistent  potentials.  They  have  chosen  Lax-Wendroff 
explicit  time  differencing.  Jensen,  and  Buot  [5]  have  used  Cayley  time  evolution  operator. 

Recent  work  has  focused  on  solving  the  time -dependent  Wigner-Poisson  equations  to 
look  for  oscillatory  solutions.  Several  groups  have  studied  dynamical  instabilities  and  undertaken 
the  simulation  the  RTD  time-dependent  behavior  using  the  discrete  Wigner  function  approach  to 
look  for  oscillatory  solutions  [5-10]. 

However  there  are  also  difficulties  in  applications  of  the  Wigner  formalism.  It  is 
important  to  recognize  that  the  Wigner  function  is  defined  over  both  real  and  momentum  spaces 
(i.e.,  which  differs  from  all  other  quantum  mechanical  formalisms)  and  this  significantly 
increases  the  theoretical  difficulty  associated  with  quantifying  the  exact  requirements  on  the 
numerical  discretization  so  as  to  insure  that  all  necessary  physical  laws  are  satisfied  (e.g.,  charge 
continuity,  momentum  balance,  detailed  balance  of  the  equilibrium  state,  and  stability  of  non¬ 
equilibrium  states).  Furthermore,  the  more  refined  numerical  solutions  of  the  Wigner-Poisson 
model  yields  results  in  disagreement  with  prior  results  that  utilized  more  coarse  grained 
discretizations. 

This  proposed  research  project  studies  the  discretization  issues  associated  with  the 
stability  and  physical  accuracy  of  the  numerical  solutions  to  the  Wigner-Poisson  modeling  of 
resonant  tunneling  diodes.  Here  the  goal  is  to  establish  the  qualitative  and  quantitative 
requirements  for  the  phase  space  discretizations  that  will  be  able  to  generate  physically  accurate 
numerical  solutions  to  the  Wigner-Poisson  model  for  open  systems.  The  goals  of  this  project  is 
also  to  study  stability  and  physical  accuracy  of  numerical  solutions  to  the  Wigner-Poisson  model. 

2.  Discretization. 

Let  us  consider  in  more  detail  a  problem  of  discretization.  An  equation  for  the  Wigner 
function  can  presented  as  follows 

V f(X,k,t)  =  T(f(X,k,t))  +  P(f(X,k,t))  +  S(f(X,k,t ))  (1) 

dt 

Here,  T  represents  the  kinetic  energy  for  electrons  in  the  device,  P  represents  the 

potential  energy  of  the  electrons,  and  S  the  effects  of  collisions  in  the  device. 

In  the  Wigner-Poisson  model  eq.(l)  has  to  be  solved  self  consistently  with  the  Poisson 
equation  to  account  for  space  charging  effect. 

The  Wigner  function,  f,  is  a  result  of  the  Weyl-  Wigner  transform  of  the  density  matrix 

/ (X,k)  =  J dsexp(iks)p(X  +^x  ")  (2) 


It  is  possible  to  write  reverse  transform 
—  jdkexp(-iks)f(X,k)  =  p(X+-,X-^) 


2  n 


(3) 


So  the  following  equation  has  to  be  valid 

f(X,k)  =  J dqf(X,q)j dsexp(is(k  -  q)) 


(4) 


If  we  introduce  discretization  a  fixed  spatial  with  a  mesh  of  we  can  rewrite  this  equation  as 
follows 


N  Ax  r  „  N  AxAk 

f(X,k)  =  —^—jdqf(X,q)Sktl  =  — — — f(X,k) 


2n 


2  n 


(5) 


As  a  result  we  have  the  relationship  between  momentum  and  spatial  intervals 


AkAx  ■■ 


71 

N  v 


(6) 


which  was  used  in  all  simulations. 

The  discrete  Wigner  function  is  periodic  in  momentum  space,  with  a  period  of 
7th  7thN  r 


Ax  L 


The  spatial  interval  is  equal  to 
Ax  = 

Nr 


(8) 


Therefore  the  following  relations  can  be  written  for  momentum  interval 
7tN .. 


A  k 


LN , 


(9) 


Hence  the  maximal  wave  vector  is  equal  to 


k 


max 


x(Nx  -1) 
2  L 


(10) 


So  it  corresponds  to  the  maximal  kinetic  energy  in  the  momentum  space. 


E 


max 


fr_ 
2 mn 


7t2h\Nx- 1)2 
8  mnL2 


(11) 


These  relations  are  taken  from  William  R.  Frensley  ]3] . 

However  the  more  refined  numerical  solutions  of  the  Wigner-Poisson  model  yields 
results  in  disagreement  with  prior  results  that  utilized  more  course  grained  discretizations. 
Simulation  runs  were  also  performed  to  show  that  the  size  of  the  momentum  mesh  can  affect  the 


shape  of  the  steady-state  current  -voltage  curve  and  the  shape  of  the  Wigner  distribution,  but  the 
size  of  the  spatial  mesh  does  not  affect  either.  The  left  plot  of  Figure  1  shows  results  using  a  fixed 
spatial  mesh  of  Nx  =  513  (where  Nx  represents  a  measure  of  spatial  grid  points  used  as  was 
defined  earlier)  and  a  variety  of  momentum  meshes,  whereas  the  right  plot  shows  a  fixed 
momentum  mesh  using  Nk  =  1024  and  various  spatial  meshes.  The  maximum  value  of  the  current 
density  and  the  shape  of  the  peaks  in  the  current  voltage  curve  change  significantly  as  the 
momentum  grid  were  refined,  but  remain  unaffected  by  changes  to  the  spatial  grid. 


Figure  1.  Steady-state  IV  curves  for  (left  side)  fixed  spatial  mesh,  Nx=5 1 3,  and  various 
momentum  meshes,  and  (right  side)  fixed  momentum  mesh,  Nk=1024,  and  various  spatial 
meshes. 


However  we  will  be  able  to  demonstrate  that  the  number  of  both  intervals  have  to  be 
equal  to  each  other  if  we  continue  to  follow  the  way  which  has  permitted  to  us  to  establish  the 
relationship  between  momentum  and  spatial  intervals. 

A  similar  derivation  can  be  made  for  the  density  matrix.  The  following  equation  has 
to  be  valid 

p{X+^,X-^-)  =  -±-\ d{p(X  + 1 , X  - 1 )\dk exp (/*(£ - 5))  (12) 

Z  Z  ZJL  Z  Z 


If  we  introduce  discretization  we  can  rewrite  this  equation  as  follows 


Mx+fx-fr 


+  i,X-t-)Sfs 
2  2  f 


As  a  result  we  have  the  second  relationship  between  momentum  and  spatial  intervals 
AkAx  =  —  (14) 

Nk 


Comparing  two  relationships  we  have  to  make  a  conclusion  that  momentum  and  spatial 
intervals  are  not  independent.  The  number  of  spatial  intervals  is  equal  to  the  number  of  intervals 
in  the  momentum  space 
Nk  =  Nx 


(15) 


In  addition  we  obtain  that  the  step  size  in  the  momentum  space  does  not  depend  on  the  number  of 
intervals 


(16) 


At  last  a  more  drastic  conclusion  has  to  be  made.  The  maximal  energy  grows  with  the 
increase  of  the  number  of  intervals.  However  all  variations  went  to  zero  before  reaching  the 
edges  of  the  wave  vector  space.  As  well  as  the  step  size  in  the  momentum  space  does  not 
depend  on  the  number  of  intervals  the  number  of  relevant  intervals  in  the  momentum  space  does 
not  also  depend  on  the  total  number  of  intervals.  So  even  in  this  case  it  is  natural  to  end 
simulation  at  lower  values  of  the  wave  vector  when  all  variations  went  to  zero  before 
reaching  the  edges  of  the  wave  vector  space  . 

In  earlier  papers  it  was  taken 

Nx  =  86  Nk  =  72  L  =  550A 


Hence  Emax=3.1  eV  is  too  large. 

The  number  of  relevant  points  in  the  momentum  space  is  possible  to  estimate  if  we  take  the 
barrier  height  0.3  eV  as  the  characteristic  value  for  the  end  of  the  variations. 

We  obtain  26  essential  intervals  for  positive  momentums.  Hence  the  number  of  relevant  points  is 
comparable  with  the  total  number  of  points.  However  the  situation  is  different  when  the  total 
number  of  points  is  large. 
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max 


x(Nx  -1) 
2  L 


(17) 


grows. 

At  the  same  time  the  number  of  relevant  points  in  the  momentum  space 


(2 mnV0)1'2  _  (2 mnV0)l,ZL 


1/2 


HAk 


fl7TN  r 


N„ 


(18) 


remains  almost  the  same  . 


3.  Stability  consideration. 


A  common  numerical  technique  for  solving  discretized  temporal  equations  is  explicit 
differencing.  However  an  explicit  scheme  is  unconditionally  unstable  .  Stability  of  the  Lax 
method  requires  that  error  in  the  discretized  equation  remains  bounded.  Fourier  analysis  of  the 
growth  of  error  in  these  schemes  has  led  to  the  Courant-Friedrichs-Lewy  stability  criterion. 


(19) 


where  v  is  the  velocity  of  the  fastest  component  of  the  solution. 

N.  C.  Kluksdahl,  A.  M.  Kriman,  and  D.  K.  Ferry  C.  Ringhofer  [4]  have  chosen  Lax-Wendroff 
explicit  time  differencing.  K.  L.  Jensen  and  F.  A.  Buot  [54]  utilize  Cayley  time  evolution 
operator.  However  Greg  Recine  does  not  use  Cayley  time  evolution  operator.  He  uses  Runge- 
Kutta  type  of  solver,  namely  the  ROCK4  subroutine.  It  is  known  that  ROCK4  is  an  explicit 
method.  Such  methods  are  unstable. 

Let  us  analyze  Stability  of  Cayley  time  evolution  operator.  It  is  enough  for  this  purpose  to 
take  into  account 


-  hk  ?f  hk 

Tf  = - -f  =  ~{3  f(x)  -  4  f(x  -  Ax)  +  fix  -  2Ax)}  (20) 

m  ox  2mAx 

for  k>0 
or 

-  hk  ?)f  hk 

Tf  = - -J-  =  — -{3  fix)  -  4  fix  +  Ax)  +  fix+2Ax)}  (21) 

m  ox  2mAx 

for  k<0. 


A  Cayley  time  evolution  operator  will  be 


f(t  +  At) 
or 


l  +  TAt/2 
1  -  TAt  /  2 


m 


(22) 


(1  -  TAt  /  2)  fit  +  At)  =  (1  +  TAt  /  2)  fit)  (23) 

According  to  von  Neuman  stability  analysis  we  are  looking  for  the  solution  as 
./"  =  4"  cxpiiqjAx)  (24) 

Hence 
for  k>0 

hk 

Tf  = - [1  -  exp(-/gAv)][3  -  exp(-/gAx)]£n  expiiqjAx)  (25) 

2m  Ax 

Substituting  (25)  into  (23)  we  obtain 
fikAt 

fn" 1  cxp(/r//Ar){  1  - [1  -  exp(-/c/Ar)][3  -  exp(-/c/Ar)] }  = 

4mAx 

hkAt 

=  %n  cxp(h//Ar){l - [1  -  cxp(-/iyAr)]  [3  -  cxpl-h/Ar)] }  (26) 

4mAr 


We  have  stable  algorithm  if  an  amplitude 


hkAt  ,  2  .  2  /  a  Ax .  ri  .  .  2.qAx.A  2hkAt  .  lrqAx.  ,qAx.r.  .  .  2.qAx.A. 
—  ^  1  ''  sirr (- — ) cos(— — )[l  +  4sin  (- — )]} 


I  #  I2  =  {1  +  4(^^)2  sin2 (^)[1  +  3 sin2 («^)]  - 

2mAx  2  2  mAx 


hkAt  ,2  _:_2,qAx^r1  ,  .^2^1^^sin2(^)cos(^x)[1  +  4sin2(^)]}-1 


)2  sin 2  (-i— — ) [1  +  3 sin 2 (-^— — )]  + 

2mAx  2  2  mAx 


{1  +  4( 

(27) 

is  less  than  1  for  q<  qth  where  the  threshold  value  of  q  is  equal  to 
qth=Jl!  Ax  (28) 


2 


Similar  calculation  for  k<0  gives 

,f..2  ,,  A, hkAt. 2  ■  2  ,  qlAx ,  r,  „  .  2,aAv  2hkAt  .  2aAx^  ,qAx.r.  .  .  2/qAx^. 
I  4  I  =  1  +  4 (— -) '  sin  -  (p—  1  +  3  sin 2  <?— )  +  — —  sin  A—  cos  A— )  1  +  4  sin 2  (A—)] 
2mAx  2  2  mAx  2  2  2 

,,  ..hkAt  .  2  .  2/qAx.r.  „  .  2/qAx 2hkAt  .  2/qAx^  aAxrt  A  .  2,qAx.^._l 

1  +  4  — — )  -  sin  -  (A— )  1  +  3  sin  (A—)] - —  sin  -  (A—)  cos  (A— )  1  +  4  sin 2  (A—)] 

2mAx  °  °  “A~  °  °  ° 


mAx 


(29) 


It  follows  from  our  solution  that  there  is  instability  at  large  q  at  any  time  step. 
It  is  interesting  that  the  threshold  value  of  q  does  not  depend  on  k. 

Let  us  compare  with  Lax  method. 


1  1  hk 

At  2  2m  Ax 


(30) 


or 


1  hkAt 

f  n+ 1  1  (  fn  |  £  n  \  1 1 ljli  (  fn  £  n  \ 


(31) 


According  to  von  Neuman  stability  analysis  we  have 


E,  =  cos(gAt)  _  ^^sin(gAx) 
mAx 


(32) 


The  magnitude  is  equal  to 


e  l2  2/  a  \  hkAt  .2  ■  2  .  »  .  .  hkAt  2  2/  a  \ 

£  1  =  cos“(q'Ar)  +  ( - )  sin-(gAr)  =  1-  ( - )  cos“(^Av) 

mAx  mAx 


(33) 


We  have  as  a  result  Courant-Friedrichs-Lewy  stability  criterion 


.  mAx 

At  < - (33) 

hk 


which  does  not  depend  on  q  but  depends  on  k. 


4. Conclusion. 


All  authors  used  the  discrete  Wigner  function  which  is  periodic  in  momentum  space.  The 
periodicity  was  proved  from  the  Fourier  transform  of  the  density  matrix.  The  inverse  Fourier 
transform  provides  us  with  the  important  proof  that  the  number  of  spatial  intervals  is  equal  to  the 
number  of  intervals  in  the  momentum  space.  In  addition  we  obtain  that  the  step  size  in  the 
momentum  space  does  not  depend  on  the  number  of  intervals.  As  a  result  the  number  of  relevant 
intervals  in  the  momentum  space  does  not  also  depend  on  the  total  number  of  intervals. 

The  analysis  of  the  stability  shows  that  Greg  Recine  has  used  an  unstable  algorithm.  . 
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